function [f, para_linear, meanutility, error_demand, diff_micro_m]...
    = gmmobj_demand_beer_cpu(para_nonlinear, para_linear_est, ...
    logshare, price_data, micro_m_data, ...
    W1, W2, mktid, x_rand, ...
    randomdraw0, rand_income, income_norm, ...
    cutoff, cutoff_craft_inc, ind_mkt_start, ind_mkt_end, n_trips, ...
    demeaned_x, demeaned_iv, brand_id2, geomkt, setup, ...
    T, iS2_brand, n_obs_per_brand, ind_obs_brand, ...
    iS2_geomkt, n_obs_per_geomkt, ind_obs_geomkt, ...
    mktsize_geomktyr, is_transformed, est_linear, x_rand2)

% PURPOSE: compute the gmm objective function value
para_nonlinear = para_nonlinear';
[n_obs, n_x_rand] = size(x_rand);

[randomdraw, para_income] = gen_randomdraw(para_nonlinear, n_x_rand,...
    randomdraw0, is_transformed);

n_x_rand2 = size(x_rand2, 2);        
mu = [x_rand, ones(n_obs, 1)] * randomdraw;        
mu = mu + (para_income*rand_income).*price_data/income_norm ...
    + para_nonlinear(n_x_rand+4)*rand_income/income_norm ...
    + (x_rand2*para_nonlinear(n_x_rand+5:n_x_rand+n_x_rand2+4, :)).*rand_income/income_norm;    

expmu = exp(mu); % exponentiated, simulated deviations from mean utilities
clearvars mu

[meanutility, is_conv] = invertmeanutility(logshare,expmu,ind_mkt_end,setup);

if is_conv
    % difference out the brand FE and geomkt FE using MAP
    demeaned_meanutility = demean_group(meanutility, brand_id2, n_obs_per_brand, ind_obs_brand, iS2_brand);

    diffy = 10;
    while diffy>1e-3
        prev_dy = demeaned_meanutility;
        demeaned_meanutility = demean_group(demeaned_meanutility, geomkt,...
            n_obs_per_geomkt, ind_obs_geomkt, iS2_geomkt);
        demeaned_meanutility = demean_group(demeaned_meanutility, brand_id2,...
            n_obs_per_brand, ind_obs_brand, iS2_brand);
        diffy = max(abs(demeaned_meanutility - prev_dy));
    end

    if est_linear
        para_linear = (demeaned_x'*demeaned_iv*W1*demeaned_iv'*demeaned_x)\(demeaned_x'*demeaned_iv*W1*demeaned_iv'*demeaned_meanutility);
    else
        para_linear = para_linear_est;
    end

    error_demand = demeaned_meanutility - demeaned_x*para_linear;
    g = demeaned_iv'*error_demand/size(error_demand, 1);
    f1 = g'*W1*g;


    [~, prob_i] = ind_shnorm(exp(meanutility), expmu, ind_mkt_end,mktid);

    % inside-good share
    n_smln = size(expmu, 2);
    n_geo = size(ind_mkt_start, 1)/T; % # of geographic mkt/year combinations

    cumsum_prob = cumsum(prob_i);
    inside_prob = diff([zeros(1, n_smln); cumsum_prob(ind_mkt_end, :)]);

    %(1) expected quantity conditional on purchasing beer
    p0 = permute(reshape(inside_prob', [n_smln, T, n_geo]), [3 1 2]); % n_geo x n_smln x T
    prob_inside_yr = 1 - prod((1-p0).^n_trips, 3);
    purch_mkt = sum(sum(p0, 3)*n_trips, 2)./sum(prob_inside_yr, 2);
    mkt_select_prob = mean(prob_inside_yr, 2).*mktsize_geomktyr;
    mean_ttq = sum(purch_mkt.*mkt_select_prob)./sum(mkt_select_prob);

    % inside good share, flavor
    mean_ttq_dummy = nan(n_x_rand, 1);
    mq_flav = nan(n_x_rand, 1);
    mi_flav = nan(n_x_rand, length(cutoff_craft_inc));

    rand_income_mkt_T_geo = rand_income(ind_mkt_start, :); % # geographic mkt/year/month combinations   x n_smln
    rand_income_mkt = rand_income_mkt_T_geo(1:T: T*(n_geo-1)+1, :);% # geographic mkt/year combinations x n_smln

    cutoff_lb = [-Inf, cutoff];
    cutoff_ub = [cutoff, Inf];
    n_cutoff = length(cutoff_lb);

    cumsum_prob_craft = cumsum(x_rand(:, 3).* prob_i);  % (craft dummy) x prob_i
    inside_prob_craft = diff([zeros(1, n_smln); cumsum_prob_craft(ind_mkt_end, :)]);

    p0_craft = permute(reshape(inside_prob_craft', [n_smln, T, n_geo]), [3 1 2]);
    prob_craft_yr_sum = sum(1-prod((1-p0_craft).^n_trips, 3), 2);
    prob_craft_yr_sum(prob_craft_yr_sum<1e-5) = 1e-5;

    purch_mkt_dummy_craft = zeros(n_geo, n_x_rand);
    for nr = 1 : n_x_rand
        cumsum_prob_nr = cumsum(x_rand(:, nr).* prob_i);
        inside_prob_nr = diff([zeros(1, n_smln); cumsum_prob_nr(ind_mkt_end, :)]);

        % (2) expected quantity_f conditional on purchasing beer
        p0_nr = permute(reshape(inside_prob_nr', [n_smln, T, n_geo]), [3 1 2]);
        purch_mkt_nr = sum(sum(p0_nr, 3)*n_trips, 2)./sum(prob_inside_yr, 2);

        %(3) expected quantity_f conditional on purchasing beer of flavor f
        prob_nr_yr = 1 - prod((1-p0_nr).^n_trips, 3);
        purch_mkt_nr_dummy = sum(sum(p0_nr, 3)*n_trips, 2)./sum(prob_nr_yr, 2);
        purch_mkt_nr_dummy(isnan(purch_mkt_nr_dummy)) = 0;
        purch_mkt_nr_dummy(isinf(purch_mkt_nr_dummy)) = 0;

        %(4) expected quantity_f conditional on purchasing craft beer
        if nr ~= 3
            for t = 1 : 12
                other_t = setdiff(1:12, t);
                purch_mkt_dummy_craft(:, nr) = purch_mkt_dummy_craft(:, nr) + ...
                    n_trips*sum(p0_nr(:, :, t).*(1 -  (1 - p0_craft(:, :, t)).^(n_trips-1) .*prod((1 - p0_craft(:, :, other_t)).^n_trips, 3) ), 2)./prob_craft_yr_sum;
            end
        end

        %(5) prob of (inc < cutoff) conditional on purchasing beer of flavor f
        income_mkt_nr = nan(size(rand_income_mkt, 1), length(cutoff_craft_inc));
        for nc = 1 : length(cutoff_craft_inc)
            income_mkt_nr(:, nc) = sum((rand_income_mkt<cutoff_craft_inc(nc)).*prob_nr_yr, 2)./sum(prob_nr_yr, 2);
        end
        income_mkt_nr(isnan(income_mkt_nr)) = 0;
        income_mkt_nr(isinf(income_mkt_nr)) = 0;

        mkt_select_prob_nr_dummy = mean(prob_nr_yr, 2).*mktsize_geomktyr;

        % aggregate across markets
        mean_ttq_dummy(nr) = sum(purch_mkt_nr.*mkt_select_prob)./sum(mkt_select_prob); %(2) expected quantity_f conditional on purchasing beer
        mq_flav(nr) = sum(purch_mkt_nr_dummy.*mkt_select_prob_nr_dummy)./sum(mkt_select_prob_nr_dummy); %(3)  expected quantity_f conditional on purchasing beer of flavor f
        mi_flav(nr, :) = sum(income_mkt_nr.*mkt_select_prob_nr_dummy)./sum(mkt_select_prob_nr_dummy); %(5)  prob of (inc < cutoff) conditional on purchasing beer of flavor f

        if nr == 3
            mkt_select_prob_craft = mkt_select_prob_nr_dummy;
        end
    end

    mq_flav_craft = sum(purch_mkt_dummy_craft.*mkt_select_prob_craft)./sum(mkt_select_prob_craft);    % (4) weighted expected quantity_f conditional on purchasing craft beer
    mq_flav_craft(3) = mq_flav(3);

    %(6) expected quantity conditional on purchasing beer and income in a bin
    %(7) expected spending conditional on purchasing beer and income in a bin
    mq_income = nan(n_cutoff, 1);
    mp_income = nan(n_cutoff, 1);
    rand_income_mkt = permute(reshape(rand_income_mkt_T_geo', [n_smln, T, n_geo]), [3 1 2]);% n_geo x n_smln x T

    cumsum_price = cumsum(prob_i.*price_data);
    avgp_mkt_unnorm = diff([zeros(1, n_smln); cumsum_price(ind_mkt_end, :)]);

    for nc = 1 : n_cutoff
        ind_rand_income_nc = (rand_income_mkt>cutoff_lb(nc) & rand_income_mkt<=cutoff_ub(nc));
        p0_nc = permute(reshape(inside_prob', [n_smln, T, n_geo]), [3 1 2]).*ind_rand_income_nc;
        prob_inside_nc_yr = 1 - prod((1-p0_nc).^n_trips,3);
        purch_mkt_nc = sum(sum(p0_nc, 3)*n_trips, 2)./sum(prob_inside_nc_yr, 2);
        mkt_select_prob_nc = mean(prob_inside_nc_yr, 2).*mktsize_geomktyr;

        purch_mkt_nc(isnan(purch_mkt_nc)) = 0;
        purch_mkt_nc(isinf(purch_mkt_nc)) = 0;
        mq_income(nc) = sum(purch_mkt_nc.*mkt_select_prob_nc)./sum(mkt_select_prob_nc);

        ind_rand2 = (rand_income_mkt_T_geo>cutoff_lb(nc) & rand_income_mkt_T_geo<=cutoff_ub(nc)); % # geographic mkt/year/month combinations   x n_smln
        price_nc = avgp_mkt_unnorm.*ind_rand2;

        p0_mkt_nc = inside_prob.*ind_rand2;
        price_mkt_nc = mean(reshape(sum(price_nc, 2)./sum(p0_mkt_nc, 2), [T, n_geo]), 1)';

        price_mkt_nc(isnan(price_mkt_nc)) = 0;
        price_mkt_nc(isinf(price_mkt_nc)) = 0;
        mp_income(nc) =  sum(price_mkt_nc.*mkt_select_prob_nc)./sum(mkt_select_prob_nc);
    end

    mq = [mean_ttq; mean_ttq_dummy];   %(1) expected quantity conditional on purchasing beer, %(2) expected quantity of type-f beer conditional on purchasing beer
    mi_flav2 = mi_flav(3, :) - mean(mi_flav([1 2 4:length(mi_flav)], :), 1);  % Pr of (inc < cutoff) conditional on purchasing craft - avg pr of (inc < cutoff) conditoinal on purchasing lager, light, foreign and ale

    mq_flav2 = mq_flav([1 2 4:length(mq_flav)]); %  expected quantity of type-f conditional on purchasing type-f where f = lager, light, foreign and ale
    mp_income2 = mp_income(2:end) - mp_income(1);  % (expected avg price conditional on purchasing beer and income in a bin) - (expected avg price conditional on purchasing beer and income in the first bin)
    mq_income2 = mq_income(2:end) - mq_income(1);  % (expected quantity conditional on purchasing beer and income in a bin) - (expected quantity conditional on purchasing beer and income in the first bin)

    % collect Micro moments:
    % •	Expected quantity conditional on purchasing beer
    % •	Expected quantity of type-f beer conditional on purchasing beer, f = lager, light, ale, foreign, craft
    % •	Expected quantity of type-f beer conditional on purchasing type-f beer, f = lager, light, ale, foreign
    % •	Expected quantity of type-f beer conditional on purchasing craft beer, f = lager, light, ale, foreign
    % •	(Expected quantity conditional on purchasing beer and income in a bin) – ((Expected quantity conditional on purchasing beer and income in the lowest bin))
    % •	[Probability of (income < cutoff) conditional on purchasing craft] – [Avg of probability of (income < cutoff) conditional on purchasing lager, light, ale, foreign]
    % •	(Expected avg price conditional on purchasing beer and income in a bin) – ((Expected avg price conditional on purchasing beer and income in the lowest bin))
    micro_m_sim = [mq(1); mq_flav2; mq_flav_craft'; mq_income2; mi_flav2(:); mp_income2];
    diff_micro_m = micro_m_sim-micro_m_data;

    if est_linear
        force_para_p = max(para_income*max(max(rand_income)))/income_norm + para_linear(1);
        f = gather((f1 + diff_micro_m'*W2*diff_micro_m)*n_obs + (force_para_p>0)*force_para_p*1e12);
    else
        f = gather((f1 + diff_micro_m'*W2*diff_micro_m)*n_obs);
    end

    fprintf(1, 'para_nonlinear='); fprintf(1, ' %1.4e', para_nonlinear); fprintf(1, ', f = %1.4e\n', f);


else
    f = 1e12;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunctions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [meanutility, is_conv] = invertmeanutility(logshare, expmu,ind_mkt_end,setup)
%% invertmeanutility
iter_count = 0; % count the number of contraction mappings

meanutility = setup.meanval_logit;
while iter_count < 100
    iter_count = iter_count + 1;
    
    exp_invert = ind_shnorm(exp(meanutility), expmu, ind_mkt_end,setup.mktid);
    
    if any(isinf(exp_invert))
        exp_invert(isinf(exp_invert)) = 1e23;
    end    
    if any(isnan(exp_invert))
        exp_invert(isnan(exp_invert)) = 1e23;
    end    
    if any(exp_invert<0)
        exp_invert(exp_invert<0) = 1e-23;
    end
    
    meanutility1 = meanutility + (logshare - log(exp_invert));
    t = abs(meanutility1 - meanutility);
    if max(t) <= setup.tol_berry_inversion; break; end
    meanutility = meanutility1;
end

is_conv = (iter_count < setup.n_iter);

function [meanprob, prob_i, exp_u] = ind_shnorm(expmeanval,expmu,ind_mkt_end,mktid)
%% IND_SHNORM
% This function computes the "individual" probabilities of choosing each brand.
% The probabilities are those associated with the normally-distributed r.c. logit.
%
% ARGUMENTS:
% expmeanval = vector of exponentiated mean utilities, sorted by market then product
% expmu = matrix of exponentiated draws of deviations from mean utilities, sorted by market then product
%
% OUTPUT:
% meanprob = vector of expected market shares, sorted by market then product


exp_u = bsxfun(@times, expmeanval, expmu);

cumsum_exp = cumsum(exp_u);
denom_exp_u = diff([zeros(1, size(expmu, 2)); cumsum_exp(ind_mkt_end, :)]);
% toc;

denom_exp_u = denom_exp_u + 1;
prob_i = exp_u./denom_exp_u(mktid,:);

prob_i(isinf(expmu)) = 1;
meanprob = mean(prob_i, 2);
